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ABSTRACT 


Reviews  and  critiques  of  selected  methods  for  nuclear  shield¬ 
ing  calculations  are  presented.  Those  techniques  considered 
are  numerical  integration  of  the  Boltzmann  equation,  moments 
method,  Monte  Carlo,  method  of  successive  scatterings,  and 
removal  cross-section  method.  An  outline  of  the  advantages 
and  disadvantages  of  each  particular  method  is  included, 
along  with  a  new  simplified  calculation  procedure. 
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1,0  INTRODUCTION 

This  report  consists  of  a  review  and  critique  of  several  selected 
methods  of  shielding  calculations .  Certain  of  the  methods  are  applicable 
to  both  neutron  and  Y"ray  penetration  problems  whereas  others  are  unique 
for  a  given  type  of  radiation.  For  each  method  outlined,  a  discussion  of 
particular  advantages  and  disadvantages  is  included. 

The  review  is  limited  to  techniques  suitable  for  shielding  studies  of 
compact  power  systems.  Certain  methods  of  analysis  are  not  considered 
because  of  their  inapplicability  to  the  problem  under  consideration.  The 
shielding  requirements  for  mobile  power  systems  such  as  aircraft  and  rockets 
place  a  premium  upon  shield  size  and  weight  in  contrast  to  stationary  nuclear 
power  plants  where  economies  is  usually  the  most  important  design  parameter. 
The  methods  of  shield  analysis  for  compact  power  systems  must  be  of  higher 
order  accuracy  than  usual  since  over  design  of  the  shield  is  undesirable. 

A  further  characteristic  of  shielding  for  mobile  systems  is  to  be  noted, 
namely,  the  geometric  configurations  encountered  will  frequently  lack  the 
symmetries  present  in  stationary  plants.  As  always,  complex  geometries 
compound  the  nature  of  the  calculational  problem  considerably. 

The  problem  of  reactor  shielding  is  concerned  with  the  biological 
shielding  of  personnel,  the  shielding  of  equipment  from  radiation  damage, 
and  the  protection  of  equipment  and  structures  from  thermal  damage.  A 
complete  discussion  of  all  the  problems  is  beyond  the  intent  of  this 
memorandum.  However,  the  fundamental  quantities  of  interest  in  any  shield¬ 
ing  analysis  are  the  neutron  and  y-ray  flux  densities  as  functions  of 
position,  energy,  and  angle.  From  these  quantities,  and  certain  physical 
parameters,  all  of  the  physical  effects  of  radiation  may  be  calculated. 

The  bulk  of  this  report  will  be  concerned  with' methods  of  computing  the 
neutron  and  Y~ray  flux  densities. 

In  the  next  section  a  brief  discussion  of  the  transport  of  neutrons  and 
y-rays  is  given.  The  third  section  of  the  report  is  devoted  to  outlining 
the  various  methods  of  use  for  shielding  calculations.  Section  4.0  contains 
a  proposal  for  a  new  method  of  performing  shielding  studies.  The  proposal 
is  untested  for  shielding  calculations.  The  last  section  is  a  bibliography 
of  papers,  etc.,  giving  a  more  general  discussion  of  shielding  theory  methods. 
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2.0  THE  TRANSPORT  EQUATION  FOR  NEUTRONS  AND  PHOTONS 
2 . I  Neutron  Transport 

The  basic  equation  of  neutron  conservation  is  the  linearized 
Boltzmann  equation.  The  neutron  flux  density,  N(r,E,G),  is  defined  as 

N(r,E,G)  =  the  number  of  neutrons  in  unit  volume  at  r; 

within  unit  energy  at  E,  within  unit  solid 
angle  about  G,  which,  in  unit  time,  cross  a 
unit  element  of  area  with  normal  in  the 
direction  G. 

The  statement  of  neutron  conservation  then  equates  losses  out  of  a  unit 
volume  of  phase  space  by  convection  and  collision  with  gains  by  collision 
of  neutrons  in  other  elements  of  phase  space  and  sources.  The  transport 
equation  is  then 


0  -  V  N(r,E,0)  +  a  (r ,E)  N(r,E,G) 


J  cr(r,E'3Q'  ;■  r,E.,fi)  NCrjE^Q') 
O' 


dfi'dE’  +  S(r,E,fi) 


(1) 


The  left-hand  side  gives  the  loss  by  convection  and  collision.  <?(r,E)  is 
the  macroscopic  total  cross-section  and  is  not  a  function  of  neutron  direction 
for  isotropic  media.  (7(r^,E,,Qf ;  r,E,G)  is  the  macroscopic  transfer  probability 
of  neutrons  from  Ef,  G1  to  E,G.  The  transfer  probability  depends  upon  the 
nature  of  the  scattering  (elastic  or  inelastic)  and  the  scattering  law.  It 
is  convenient  to  consider  the  Boltzmann  equation  in  terms  of  the  lethargy 
variable  u  =  ln(EQ/E),  where  E0  is  some  peak  energy.  For  fission  sources., 

E0  is  usually  taken  to  be  15-20  Mev.  For  later  illustrations,  the  transport 
equation  for  slabs  will  be  used.  For  slab  geometry,  and  in  terms  of  the 
lethargy,  Equation  (1)  is  written 


^N(x,u,!i)  +  o  (x,u)N(x,u,|J,) 


J  J1 


a  (x,u 1  ,}Jt 1 ;  x,u,}Ji)N(x,u  1  ?\i 1 )  du1  d^ 1 


+  S(x,u,jj0 


(2) 


where  ^  =  cos  9,  with  0  the  angle  between  the  +x  axis  and  the  direction  of 
neutron  motion. 
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Various  special  scattering  laws  are  easily  treated.  If  all 
scattering  is  assumed  spherically  symmetric  in  the  laboratory  system, 
then  the  transfer  probability  is  independent  of  and  |A.  To  first  order 
accuracy  inelastic  scattering  is  spherically  symmetric  in  laboratory 
coordinates.  For  elastic  scattering  symmetric  in  the  center  of  mass 
system,  the  lethargy  and  scattering  angle  are  uniquely  related.  Hence 
the  integral  over  u1  and  [X1  involves  a  6  function  in  one  or  the  other 
variable.  Note  that  even  if  the  transfer  probability  is  independent  of 
1  ,  p,  the  flux  density  is  still  a  function  of  angle. 

Boundary  conditions  for  the  transport  equation  for  shielding 
problems  are  usually  given  in  terms  of  incident  neutrons  at  interfaces. 

Thus,  for  a  slab  shield  of  thickness  a,  adjoining  a  core  at  x  *  0,  the 
boundary  conditions  would  be  given  at  x  *  0  for  0  £  p>  &  1.  At  x  *  a, 
assuming  a  vacuum  boundary,  the  condition  is  N(a,u,ji)  *  0^  -1  ^  |i»  £  0. 

Except  for  rare  cases,  the  external  source  is  zero,  or  a  delta  function 
at  x  *  0  that  is  given  by  the  boundary  conditions . 

2 .2  y-ray  Transport 

The  transport  equation  for  photons  is  the  same  as  (1) 5  however, 
simplifications t are  possible.  The  only  source  of  scattered  photons  is  the 
Compton  effect'  For  Compton  scattering  the  energy  and  angle  are  related, 
and  again  a  delta  function  will  appear  in  the  integfand.  It  is  convenient 
to  change  variable  from  the  flux  density  to  the  energy  density,  I,  where 

I(r,E,Q)  =  E  *  N(r,E,Q) 

In  terms  of  I,  equation  (1)  becomes 

Q  *  Vi(r,E,Q)  +  CT(r,E)  I(r,E,Q)  -  J 

E' 

+  S'(r,E,£2)  (3) 

If  E  is  expressed  in  electron  mass  units,  then  the  photon  energy  may  be 
written  in  terms  of  the  Compton  wave  length  X9  with 

X  »  1/E  (E  in  electron  masses) 


j  a  CE.jEj'Q1  ;  x,E,n)  I(r,E',n')  |r  dE'dfi’ 
fi' 


(*) 


Other  sources  such  as  annihilation  radiation  and  bremsstrahlung  are 
usually  negligible. 
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with  the  above  change  of  variables,  the  scattering  kernel  is  written 
a(r,E',n';  r,E,Q)  =  a(r,X'  ,0-0')  6(1  +  A*. 

The  scattering  kernel  k(Xf,X)  is  defined  as 

k(X',X)  =  2n  p-CT(\',fi'.Q) 

and  is  given  by  the  Klein-Nishina  formula  for  Compton  scattering.  In 
particular , 

k(X',X)  =|  jp  [£t+  j^-+  2(X'  -  X)  +  <X'  -  X)2]  *  A  (X'  £  X.  s:  X1  +  2) 

*  0  otherwise,  where  A  is  the  atom  density. 

For  slabs  the  transport  equation  is  thus 

X 

V-  ^  I(*,X,|i)  +  a(x,X)  I(K,X,n)  =  J  J  I(X,X',H») 

o  n» 

+  S(l  +  X1  -X  -0*0')  dOfdXf+  S 1  (x,X  ,|i)  (4) 


The  transport  equation  for  photons  is  somewhat  simpler  than  for 
neutrons  because  of  the  known  scattering  law.  In  particular,  only  photons 
in  the  wavelength  range  X  -  2  ^  Xf  £  X  can  contribute  to  the  scattering 
source.  This  simplification  of  the  scattering  source  is  particularly 
useful  in  numerical  integration  and  moments  expansions  of  the  transport 
equation  for  photons. 
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3.0  SELECTED  METHODS  OF  SHIELD  CALCULATIONS 


In  view  of  the  complex  nature  of  the  Boltzmann  equation  analytic 
solutions  are  possible  only  for  the  simplest  of  problems.  A  number  of 
procedures  for  approximate  shield  calculations  have  been  proposed.  All 
of  the  methods  represent  approximations  to  the  transport  equation,  either 
numerical  approximations,  or  simplifying  physical  approximations.  Several 
methods  of  both  types  will  be  considered. 

3 . 1  Numerical  Integration  of  the  Boltzmann  Equation 

There  are  numerous  numerical  approximations  used  for  transport 
calculations.  Rather  than  outline  a  particular  method,  such  as  the  S^ 
method,  a  general  procedure  is  considered.  The  approach  is  sufficient  to 
indicate  the  nature  of  the  approximations  inherent  in  the  method. 


Basically  every  numerical  approximation  is  found  by  dividing 
continuous  variables  into  discrete  variables  and  solving  the  resultant 
difference  equations.  Values  intermediate  to  the  discrete  points  are 
found  by  interpolation.  For  the  transport  problem,  for  neutrons,  the 
lethargy,  angle,  and  spatial  ranges  are  discretized.  First  the  lethargy 
interval  0  £  U  ^  Ut^  is  divided  into  G  groups,  not  necessarily  of  equal 
spacing,  with  AU  =  A  =  U  -  U  ^ .  A  thermal  group  with  index  G-fl  is 
appended.  g  g  g  g  Integration  of  the  Boltzmann  equation 

(2)  over  the  interval  A  yields 


M*  ~  NS(x,p)  +  ag(x,p)  Ng(x,p0  -  S  du  ^  a(x,u 1  ;  x,u,p)  N(x,u,;p’) 


5-1 


u*  -1 


du 


SS(x5|i) 


(5) 


The  quantities  Ng,  ag,  Sg  denote  average  values  of  the  variable  in  the 

interval  A  .  Thus, 
g 


NS(x,p,) 


5-1 


N(x,u',p)  du 1 


and 


Pu 

g  a(x,u  1 )  N(x,u?,p)  du' 


ag(x,M0  = 


s  N(x,u 1  jP>)  du1 
*3 


g-i 
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Note  that  although  a(x,u)  is  not  a  function  of  p,  the  average  value  of  the 
cross-section,  ~"g,  is  a  function  of  p,.  Usually  the  variation  of  ”g  with 
is  ignored  siSce  N(x,u,p,)  is  a  reasonably  flat  function  of  ^  for  reactor 
calculations.  In  shield  calculations  there  may  be  a  very  highly  anisotropic 
flux  distribution,  particularly  at  deep  penetrations  and  high  energies ;  in 
which  case  the  distinction  is  important. 


An  immediate  difficulty  with  the  numerical  procedure  is  apparent. 

The  proper  cross-sections  for  the  group  equations  cannot  be  computed  until 
the  solution  is  known.  In  reactor  calculations  the  properly  weighted  cross- 
sections  are  frequently  found  by  assuming  a  spectrum  initially  or  by  computing 
on  approximate  spectrums  from  an  infinite  medium  calculation.  The  infinite 
medium  method  is  useful  for  shield  calculations  where  the  spectrum  nhardeningM 
with  increasing  depth  may  be  explicitly  included. 


Assuming  that  properly  weighted  coefficients  may  be  found  then  the 
scattering  integral  can  be  handled.  The  integration  over  u1  is  approximated 
by  a  sum,  such  as 


X,u,|l)  N(x,u  1  ,(J. ' ) 


Gf  1  _ _  _ 

^  cts  ,s(xj|J.';x,|j,)Ns  (x_,p,’) 
g’=l 

(6) 


where  g!,g  represents  the  transfer  probability  from  group  g!  to  group'  g, 

The  determination  of,  ~g f’,  g  depends  upon  the  scattering  law  and  the  weighting 
function  used  for  the  lethargy  integration.  The  sum  includes  all  of  the 
groups.  For  shielding  purposes  the  thermal  region  is  rarely  divided  Into 
many  groups,  hence  shield  calculations  characteristically  contain  fewer 
groups  than  criticality  calculations.  In  general,  low  energy  neutrons  are 
rapidly  absorbed  and  do  not  up-scatter.  In  this  case,  all  elements  of  g 1 ? g 
are  zero  for  g1  >  g.  For  elastic  scattering  only  elements  for  g'  near  g  are 
non-zero.  For  inelastic  scattering  g’  may  be  far  removed  from  g. 

The  angular  variation  is  treated  in  similar  fashion.  The  range 
-1  ^  |i»  ^  1  is  divided  into  segments  of  width  Ap^  =  A^  =  l  *  vari°us 

terms  in  the  group  transport  equation  (3)  may  be  integrated  over  A  as  before. 
The  integration  yields  n 

f”  “  df 1,8 ■ <*•»  *  ■£  +  i.f«] 


aS(x,pO  NS(x,ii,)  d|i>  = 


Ng(x)  +  N  8(x) 

L  n  n-1  J 


1 
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where  a  ,  |3  are  weight  coefficients  which  depend  upon  the  quadrature  rule 
assumed1}  For  instance,  in  the  S  method  TTg  .  v  is  assumed  to  vary  linearly 
over  an  interval.  For  the  Wick  method,  the 'gdussian  weight  functions  are 
used.  Obviously  many  different  possibilities  exist. 

The  transfer  probabilities  from  equation  (6)  may  be  written 


u  ^  _  _  _ _ 

J  dp,  J  dp'[^CTS  ’S  (x  ,|J. ' ;  x,p)  Ns  (x ,  p, ' )  dp,  *  =  ^  ^CTn’’n('X')  Nn'  ^ 

^n-l  **1  g1  n‘  g' 


The  group  transport  equation  is  then  of  the  form 


srk8« +  V.g«]+  e„ks« +  -I  V8'« 


n*  g’ 


+  «n«M 


(8) 


The  spatial  variation  is  treated  as  before.  In  general  the  result  is  a 

difference  equation  coupling  various  mesh  points  in  the  indices  g,n,j 

(j  is  a  spatial  index)  .  A  typical  difference  equation  would  be  of  the  form 


r 

n 


3  n,j 


-1-  1 


g 


.  .  N  .  - 

n,j  n,j-l 


+  t 


g 


.  N  1  , 
n,j  n-l,j 


+  b 


n,j 


N  .  S  = 

n-l> j-1 


I 


n 


T  f  .S? ,g  N  f  . 
n1  ,n; J  n  , j 


+  S  . 
n,j 


g 


g  -  1,2,  .  .  OH 

n  ~  0, 1 ,  .  ♦  . ,  N 
j  =  0 , 1 ,  .  .  • ,  J 


(9) 


The  result  is  a  very  large  set  of  simultaneous  algebraic  equations. 

The  method  of  solution  is  straight  forward.  An  initial  distribution, 
consistent  with  the  boundary  conditions,  is  assumed.  The  difference  relation 
(9)  is  applied  repetitively  to  each  point  until  a  convergent  solution  is 
obtained.  Further  comments  on  the  solution  of  algebraic  equations  are  given 
in  Section  4.0. 
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From  the  elementary  sketch  above,  it  is  possible  to  assess  the 
advantages  and  disadvantages  of  the  direct  integration  method.  The 
advantages  are  several: 

1)  Inhomogeneous  media  and  finite  media  are.  readily 
incorporated  into  the  analysis. 

2)  Complex  scattering  kernels  can  be  approximated  by 

the  transfer  matrix  g',g. 

cr  i 
XL  ,n 

3)  The  truncation  error  associated  with  discretizing 
can  be  reduced  by  refinement  of  the  mesh. 

4)  The  method  yields  all  of  the  desired  information, 
i  *e  . ,  N(r;u,nO  . 

5)  The  results  are  not  subject  to  statistical  variation, 
hence  the  method  is  readily  adapted  for  deep  penetrations . 

6)  Surveys  may  be  made  rapidly  by  reducing  the  number  of 
discrete  variables. 

The  disadvantages  associated  with  the  method  are  rather  severe. 

1)  For  detailed  studies  the  amount  of  machine  time  is 
very  large.  In  particular,  the  number  of  iterations 
to  reach  a  solution  is  roughly  proportional  to  the 
square  of  the  unknowns.  To  double  the  mesh  size  in 
each  dimension  increases  the  computation  time  by  a 
factor  of  64  approximately. 

2)  Higher  dimension  problems  .,  and  irregular  geometries  , 
increase  the  computational  burden.  Two  and  three 
space  dimension  problems  are  at  the  limit  of  present 
machine  capabilities. 

3)  The  proper  weighted  cross-sections  are  difficult  to 
find  and  require  a  subsidiary  calculation  of  the  same 
order  of  magnitude  as  the  original  problem. 

4)  Reduction  of  truncation  error  increases  the  computational 
time  in  roughly  quadratic  fashion. 

3.2  The  Moments  Method 


The  basis  of  the  moments  method  is  an  expansion  of  the  angular 
and  spatial  dependence  of  the  flux  density.  The  objective  of  the  expansion 
is  to  decouple  the  complex  behavior  of  the  flux  density  into  a  set  of 
variables  which  may  be  recombined  to  approximate  the  flux  density.  The 
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expansion  should  be  carried  out  in  such  a  manner  that  the  resulting 
simultaneous  equations  for  the  expansion  variables  are  easy  to  solve,  and 
further,  the  expansion  variables  should  form  a  rapidly  convergent  sequence. 

The  process  of  reduction  is  outlined  for  a  plane  problem  for 
photons .  The  transport  equation  is 


dx  +  CT(X.’^)  =  J  J  I(x,X  1  ,p/)  k  (X  1 ,  X )  * 

0  O’ 


*  6  (i  +  x  ’  -  x  -  o  •  n 1 ) 

2tt 


dQJ  dA. 1  4-  S  (x  ,p»)  . 


(10) 


The  expansion  for  the  angular  coordinates  is  in  terms  of  the 
Legendre  polynomials  Py  (p)  whereas  the  spatial  expansion  is  in  terms  of 
the  spatial  moments  moments  xn .  The  generalized  moment,  say  bn  p(A.), 
is  determined  as 


0_  nil 
2rrg0 
n! 


r>”  A 

j  dx  J  dp-P^(M-)xnI(x,X,tJ,) 

_co  _i 


(11) 


where  an  is  the  cross-section  at  the  largest  source  energy.  The  choice 
of  constants  facilitates  normalization. 

The  objective  is  to  reduce  the  transport  equation  to  coupled 
equations  (through  n  and  (?)  for  the  variable  b^  g(\). 

The  first  step  in  the  reduction  of  the  transport  equation  is  to 
carry  out  the  angular  expansion.  The  expansion 

I(x,X,n,)  =  £  ^i_±_ A  i  (x,X)  P^(ii) 

H 

is  used.  By  the  orthogonality  of  the  P j^(M') 

I^(x,X)  =  2tt  J  I(x,X,|a) 

-1 


By  multiplying  equation  (10)  by  P^ (p)  and  integrating  over  all  Q  we  have 
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V 


5li 

dx 


+ 


a(x,\)Ig  = 


j  Pj(l  +  X'  -  X)  k(X',X>  * 
0 


I^(x,X*)  dX*  +  Sjj(x,X)  (12) 

The  proof  of  (12)  follows  from  the  addition  theorem  for  Legendre  polynomials 
and  the  known  recursion  relations  for  the  P Jji) „  Equation  (12)  applies  for 

i  =  0,1,2,  .  .  .  .  k 

The  spatial  dependence  is  eliminated  from  the  definition  (11) , 

that  is 


n-t-  1 


b»,iw 


OD 

Io(x,X)  xndx 


To  find  the  equation  for  the  b  (X)  we  multiply  (12)  by  — -p —  x11  and 
integrate  over  all  x.  The  '*  left-hand  side  is  *  considered 

first.  The  derivative  terms  are  integrated  by  parts  to  yield 


rrf  1 


ft+l  go _  f  nT 

2jFI  nl  l  x  1+1 


oo 


n-1- 


1  V 


nb  1 


nx  I^.L(x,X.)dx  j  +  2j^I  lX  \-l  [ 


rtf  1 


-1  1  CT  p® 

tn‘  I._1(x,X)dx|  +  -^T —  (  a  (x,\)I. (x,X)dx 


(13) 


In  order  for  the  integrated  terms  to  vanish,  it  is  necessary  that  I*  vanish 
faster  than  x  for  large  x0  An  immediate  consequence  of  this  restriction  :i.s 
that  the  medium  must  be  infinitely  thick,  Obviously  this  represents  a 
serious  restriction.,  Equally  important  is  that  fact  that  a(x,X)  must  actually 
be  independent  of  x,  and  thus  the  moments  method  only  applies  to  homogeneous 
media.  The  above  two  restrictions  limit  the  applicability  of  the  method 
particularly  to  compact  power  systems. 

With  the  above  assumptions,  equation  (13)  becomes 


(14) 


bn>Jd)  -  o0  v.,nw  4tVi,h«} 
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The  right-hand  side  of  equation  (12)  is  easily  found  to  be 


X 

P^U  +  X'  -  X)  k(X',X)  bn  |(X!)dX!  +  Sn  j(X) 


The  Boltzmann  equation  reduces  to 


ct(X)  bn  ^(X)  =  J  P^(l  +  X.1  -  X.)  k(X',X)  b^  ^(X')dX'  + 

CTo  ilf+T  bn-l,JJfl^')  +  2jT+r  +  Sn,JL^ 


(15) 


The  aim  of  the  method  has  been  achieved  in  that  a  coupled  set  of 
simple  equations  have  been  found  for  the  moments  of  the  expansion.  The 
integral  is  evaluated  numerically  since  the  functions  P  (X  1  )  and  k(\',X) 
are  known.  Further,  the  integration  is  over  the  range  X.  -  2  £  X1  ^  X  and 
hence  the  entire  integrand  is  known. 

For  the  case  of  neutron  penetration,  the.  integration  is  more 
involved  since  the  range  of  integration  is  from  0  ^  u1,  and  furthermore, 

the  scattering  kernel  is  not  as  simple  a  function  as  for  y-rays .  For  neutrons 
an  approximate  kernel  is  derived  in  a  manner  similar  to  the  methods  of  direct 
numerical  integration. 

The  equation  for  the  moments,  (15),  is  solved  sequentially  for 
n  =  1,  2  .  .  .  .  Note  that  the  equation  couples  only  lower  moments  together 
and  hence  no  truncation  is  necessary  to  solve  for  any  given  moment.  This  is 
in  contradistinction  to  the  usual  spherical  harmonics  method  of  neutron 
transport. 


The  reconstruction  of  the  energy  density  is  reasonably  straight 
forward.  It  is  sufficient  to  find  the  coefficients  Ij^(x^X) .  The  objective 
is  to  find  a  rapidly  convergent  expansion  for  the  function  I^(x,X)'.  The 
coefficients  b  «(X)  are  not  directly  useful  since  the  factors  xl  do  not' 
obey  a  simple  3  *  orthogonality  rule.  The  following  physical  approximation 
is  usually  applied.  The  behavior  of  the  flux  density  with  increasing 
penetration  is  roughly  exponential.  An  expansion  of  the  form 


Ij^(XjX) 


I 


-CTnX 

e  °  P, 


(x) 


(16) 


with  p  (x)  a  polynomial  of  degree  n,  is  assumed.  The  coefficients  can  be 
evaluated  as 


TM-.3522 


12 


3n>^(X)  “  L  Ii(X^)PnW 


(1?) 


if  the  p  (x)  are  the  Laguerre  polynomials.  But 
b  £(X)  ,  fehe  an  j j(X)  are  linear  combinations  of 
n '  3 -  1,2,  i,  .  n.  That  is, 


from  the  definition  of 
the  b  ,  a(X)  for 

ft  ,  A, 


/ 

/ 


c  ,b 
n  n 


(18) 


Several  points  are  worth  noting.  First ,  the  expansion  (16)  was 
particularly  chosen  so  that'  the  coefficients  were  easily  found,  as  by 
equation  (17) .  Obviously  other  weight  functions  than  e“aox  are  possible. 
Different  weight  functions  give  different  coefficients,  but  orthogonality 
relations  relative  to  a  weight  function  are  readily  found.  The  choice  of 
e“aox  is  particularly  useful  since  only  a  few  a  AN  4  - 

the  series  converges  rapidly. 

The  Boltzmann  equation  could  have  been  expanded  directly  in  terms 
of  the  an  £(\)  instead  of  the  bn  £(\)  .  In  a  certain  sense  the  bnj£(X)  are 
more  general  in  that  the  appropriate  coefficients  for  weight  functions 
other  than  unity  are  most  readily  obtained  from  the  bn  £(4)  *  Thus  the 
bn  jj(^)are  the  simplest  to  use  for  finding  expressions  such  as  (18). 

The  basic  elements  of  the  moments  method  are  now  clear.  The 
calculational  procedure  is  simple  (particularly  for  photons)  and  the 
reconstruction  of  the  flux  density  allows  wide  latitude  for  use  with 
rather  involved  weight  functions.  The  development  does  illustrate  the 
particular  advantages  and  disadvantages  of  the  method . 

Advantageous  properties  of  the  method  include 

1)  The  method  yields  quite  accurate  results  for 
deep  penetrations . 

2)  Relatively  universal  penetration  curves  can  be 
derived  and  used  to  study  composite  media  with¬ 
out  additional  calculation, 

3)  Machine  time  is  relatively  small  compared  to 
direct  integration. 

4)  The  method  lends  itself  readily  to  deriving 
elementary  equations  which  approximate  the 
solution . 

3)  Ancillary  information  such  as  build-up  factors 
are  easily  found. 


\f*  J  a.  A.  C  UCCUCU  ,  JL  «  C  .  , 
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The  disadvantages  are  quite  serious  for  compact;  shields. 


1) 

2) 

3) 


4) 

5) 


The  requirement  of  an  infinite  medium.  By  uke 
of  approximate  formulas  this  restriction  can  b\ 
made  less  severe. 

Homogeneous  media.  This  is  very  serious  and 


\ 


not  removable.. 


The  determination  of  the  angular  distribution 
requires  many  approximations  since  the  proper 
weight  functions  in  the  expansions  are  not 
readily  apparent. 

The  results  are  generally  applicable  only  beyond 
several  mean  free,  paths,  thus  the  full  flux  density 
is  not  found. 

In  many  cases  the  accuracy  of  the  result  is  hard 
to  determine. 


3.3  The  Monte  Carlo  Method 


The  Monte  Carlo  method  is  a  statistical  sampling  procedure  for 
conducting  theoretical  experiments  on  particle  distributions.  In  its 
simplest  form  the  method  consists  of  a.  straightforward  calculation  of  many 
particle  histories.  The  results  yield  a  probabilistic  measure  of  the 
penetration  and  distribution  of  particles.  Usually  a  straightforward 
reproduction  of  histories  is  inefficient  and  a  variety  of  modifications 
to  the  procedure  are  adopted  to  yield  equivalent  results  with  less  effort. 

The  basis  of  Monte  Carlo  calculations  is  the  Central  Limit  Theorem 
of  Statistics..  The  theorem  relates  theoretical  first  and  second  moments  of 
distribution  functions  to  the  normal  distribution.  Let  f  (x)  denote  a 
distribution  function,  hence 

J  f(x)dx  *=  1 


The  first  moment  of  the  distribution,  the  mean  value,  is  then 


<x> 


xf  (x)  dx 


The  second  moment  is 


<x2> 


x^f  (x)dx 
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A  very  useful  statistical  quantity  is  the  second  moment  about  the  mean, 
usually  called  the  variance  r „  Evidently 


CD 

r  =  <(x  -  <x>)2>  *=  (x2  -  2x<x>  +  <x>2)  f  (x)  dx 

u 


2  2 
=  < X  >  -  <x> 


The  square  root  of  the  variance  is  called  the  standard  deviation  and 
denoted 


a  *=  /r 


The  Central  Limit  Theorem  then  states  that  if.  a  finite  number 
uf  values  of  x,  say  x^  ,  1  £  i  ^  N, 

-  1  V 
x  =  i4xi 

i 

x  -  <x> 

Then  the  variable  a  has  a  distribution  that  approaches  a  normal 

distribution  for  large  N.  In  other  words 


x  «  <x> 

with  an  approximately  normal  distribution  with  standard  deviation 
—  a 

ct  =  7n 

The  approach  to  the  normal  distribution  is  accurate  to  terms  0  (1//N)  which 
is  negligible  for  large  N. 

The  importance  of  the  Central  Limit  Theorem  is  apparent.,  since 
it  is  then  possible  to  evaluate  the  mean  of  a  distribution  and  to  find  the 
limits  of  confidence  on  the  calculated  value.  This  is  precisely  the 
problem  the  Monte  Carlo  is  addressed  to. 

Generalizations  of  the  Central  Limit  Theorem  are  readily  found „ 
Thus,  if  g(x)  is  some  function  of  x,  where  x  is  distributed  according  to 
the  distribution  f (x)  then 
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and 


r°° 

<g>  =  g(x)f  (x)dx 


„  p”  ~ 

<g  ?  =  g  (x)f(x)dx 

d 


and 


rg  =  <%>  -  <g>2 


The  Central  Limit  Theorem  is  then 


g  -  ££  g(x.)  W  <g> 


g  /K 


<KV  -  <£>' 


In  most  cases  the  Monte  Carlo  method  is  used  to  find  g  with  a 
not  known.  A  statistical  estimate  of  a  may  be  found  readily.  The 
expected  value  of  g  is  6 


g 


CO 

<s>  ’III  g(x.)f(x)dx  =  |  U  <g>]  =  <g> 


that  is  the  expected  value  of  g  is  the  mean.  Hence  g  is  an  unbiased 

estimate  of  <g>  .  The  expected  value  of  the  sample  variance  ~~2  —2  is 

g  -  g 


<g2-g2>“j|J  ^  (g(x..)  -  g)2  f (x) dx 

i 


1 

N 


£  J  f  (x) dx  |^(g(xi)  -  <g>)2 - 


(g  -  <g»2] 
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2 

a 

g 


n  r 

N-l  L8 


2 


-2 

-  g 


Thus , 


Hence  a  statistical  estimate  of  the  standard  deviation  for  the  distribution 
of  g  is  readily  available. 

The  two  central  problems  of  the  Monte  Carlo  method  are  seen  to 
be : 


1)  Selection  of  x,.  from  a  distribution  f(x) 

2)  Reduction  of  the  sample  variance 

The  selection  from  a  distribution  function  is  usually  accomplished 
by  use  of  random  numbers,  i<,e.,  numbers  distributed  uniformly  in  the  interval 
0  £  §  £  !•  If  the  cumulative  distribution  F(x)  is  defined  as 


f  (x)  dx 


and  if  a  random  number  §  is  chosen  then 
•  X  =  F_1(5) 


is  a  random  variable  x  selected  from  f  (x) .  This  follows  since 


5  =  F(x) 


and 


P  (§)d§  =  p  (x)  dx 
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Since  p(§)  *  1  we  have 

p(x)  ■  ^  “  f(x) 

hence  x  distributed  according  to  f (x) 

For  multi -dimensional  and  complicated  distributions  direct 
evaluation  may  be  impractical  and  the  rejection  technique  is  often  used. 
The  method  is  illustrated  for  a  one  dimensional  distribution  f  (x)  .  If 
f (x)  has  the  distribution  as  shown  below 


where  the  limits  are  a  £  x  ^  b  and  f  =  max  f  (x)  ,  then  the  rectangle  abed 
contains  all  of  f  (x)  .  Two  random  numbers  §  ,  ^2  are  used  to  determine  a 
point  in  the  rectangle.  If  (5p  ^es  f  00  then  is  accepted, 

otherwise  the  process  is  repeated*  The  geometric  interpretation  of  the 
method  Is  evident , 

There  are  many  other  methods  of  selection  which  are  not 
enumerated  here . 

The  reduction  of  sample  variance  may  obviously  be  achieved  by 
increasing  the  number  of  samples  N,  Alternate  methods  are  to  modify  the 
sampling  probabilities  to  decrease  the  sample  variance  2  —2*.  Particularly 

u-seful  techniques  are  considered  for  the  penetration  ®  ®  problem. 

Let  p  denote  the  probability  of  a  particle  penetrating  the  shield 
and  (1-p)  the  probability  of  not  penetrating.  If  N  particles  are  incident, 
then  the  probability  of  S  penetrating  is 

p  Nl  Sn  nN-S 

Ps  “  Sl(N-S)  :  p  U-P) 


i.e.,  a  binomial  distribution.  The  mean  of  the  distribution  is  Np  and 
the  variance  Np(l-p).  The  fraction  transmitted  is  <s>  with  variance 

p  (1-p)  The  error  in  the  estimated  value  of  p,  N  say  e,  is 

N 
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and  the  fractional  error  is 


(19) 


Obviously  for  small  p  the  fractional  error  will  be  small  only  for  large  No 

To  reduce  the  error,  i.e.,  the  sample  variance,  it  is  frequently 
possible  to  modify  the  penetration  probability  in  a  known  manner  such  that 
the  new  penetration  probability,  say  p',  is.  greater  than  p.  In  particular, 
if  p1  can  be  made  equal  to  unity  then  the  error  is  zero  and  p  can  be 
determined  exactly.  Note  that  the  transformation  from  p  to  p '  must  be 
Known  so  that  p  may  be  computed  from  the  estimated  p1.  An  alternative 
procedure  is  to  modify  the  calculation  procedure  such  that  the  correct  p 
is  estimated  but  the  sample  variance  of  the  procedure  is  different  from 
(19),  and  hopefully  smaller.  Examples  of  both  types  of  reduction  of 
variance  will  be  shown. 


The  increase  of  penetration  probability  may  be  achieved  by  the 
exponential  transformation.  The  transformation  is  best  understood  by 
reference  to  the  Boltzmann  equation  in  the  form 


>*  §+ oi 


F(x,X,M-> 


The  spatial  variation  of  I(x)  is  of  the  form 
~x 

I(x)  =  e  ^  jjs  (x,X  ,|i) J 


For  straight  ahead  particles,  i.e.,  jjl  =  1,  the  penetration  is  roughly 
proportional  of  e"0^.  For  thick  shields  the  probability  of  penetration 
is  very  small.  To  increase  the  penetration  probability  the  transformed 
variable 


-  (Vy 

I(x,X,H)  =  e  I(x,\,|A) 


is  considered.  In  terms  of  the  transformed  variable  the  Boltzmann  equation 
is 
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V>  1  +  (.a  -  a)  l  =  F'(x,X,M<) 

* 


and  hence  the  attenuation  factor  may  be  reduced  by  proper  choice  of  O'. 
Usually  Qt  is  slightly  smaller  than  the  minimum  of  cr.  The  Monte  Carlo 
evaluation  of  I  yields  an  estimate  p'  with  small  fractional  error, 
p  is  then  merely  e"^xp’  with  the  same  fractional  error. 

Other  variance  reduction  schemes  can  be  used  with  yield  the 
same  p,  but  different  variance.  A  particularly  notable  method  is  the 
importance -samp ling  technique.  As  applied  to  penetration  problems  the 
basis  of  the  method  is  the  following.  When  a  particle  is  heading  out  of 
the  shield  it  is  "important"  to  follow  the  particle  since  it  may  penetrate. 
Conversely,  a  particle  heading  back  into  the  shield  has  a  smaller  chance 
of  penetrating  and  hence  is  less  important.  It  is  advantageous  to  bias 
the  particle  behavior  so  that  collisions  favor  forward  scattering  of  the 
particle  at  the  expense  of  back-scattering.  Thus  the  probability  of 
forward  scatter  might  be  made  n  times  more  likely  than  back-scatter.  To 
eliminate  the  bias  introduced,  the  forward  scattered  particle  is  not 
counted  as  one  particle  but  1/nth  of  a  particle.  The  result  of  the 
sampling  procedure  is  that  most  of  the  computation  time  is  spent  treating 
forward  moving  particles. 


Extensions  of  the  method  are  readily  available.  Since  the  bias 
scattering  introduces  a  ’‘weight"  associated  with  particles  we  may  use  other 
weighting  procedures.  For  instance,  instead  of  terminating  a  particle 
history  by  absorption,  we  might  reduce  the  particle  weight  by  c 7a/at  at 
each  collision  and  make  all  collisions  scatterings.  In  this  case,  no 
particle  dies  by  collision  but  only  when  the  weight  gets  too  small  to 
contribute  sensibly  to  the  penetration. 


A  further  extension  of  the  same  basic  procedure  is  the  "Splitting 
and  Russian  Roulette"  procedure.  In  this  method  one  emphasizes  deep  pene¬ 
tration  by  splitting  each  particle  that  penetrates  to  a  given  depth  into 
several  particles  of  reduced  weight.  If  a*particle  returns  toward  the 
source  then  a  subsidiary  game  of  Russian  roulette  is  played,  where  the 
particles  chance  of  survival  is  equal  to  the  reduced  weight  of  the  splitting. 
Usually  several  splitting  planes  are  introduced  at  various  depths. 


Another  process  of  using  weights  is  the  use  of  expected  values. 
Suppose  a  particle  is  at  point  x  heading  in  the  direction  of  .  The 
probability  of  penetration,  without  further  collision  is  calculable  as 


W  = 
o 


--(T-x) 


with  T  the  shield  thickness.  An  amount  of  penetration  Wq  is  tallied  and 
the  particle  is  given  the  reduced  weight  (1  -  W  ) .  The  particle  is  then 
followed  until  a  collision  occurs  or  the  particle  penetrates.  If  the 
particle  penetrated  the  score  (1  -  WQ)  is  tallied  and  a  new  particle 
studied.  If  a  collision  occurs  the  new  energy  and  direction  parameters 


t 
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are  determined  (and  perhaps  a  new  weight) .  Let  the  resulting  particle 
have  a  weight  0^(1  -  WQ) .  The  penetration  without  further  interaction 
is  W^.  The  score  0^  (1  -  WQ)  is  tallied  and  the  remaining  particle 
is  given  the  weight  Q^(l  -  WQ) (1  -  W^) .  The  process  continues  until  the 
weight  becomes  negligible  or  the  particle  penetrates. 

Occasionally  a  particle  with  very  small  weight  reaches  a  very 
important  region  of  the  problem  -  for  instance  very  near  the  edge  of  the 
shield.  In  such  a  case,  the  particle  need  not  be  killed  because  of  the 
low  weight,  but  may  be  given  a  chance  to  survive  by  the  Russian  roulette 
game.  If  the  particle  survives  the  Russian  roulette,  the  weight  is 
increased  accordingly.  This  procedure  is  obviously  a  combination  of 
importance  sampling  and  Russian  roulette. 

The  estimated  statistical  variance  when  using  any  or  all  of  the 
above  methods  is  frequently  difficult  to  ascertain.  For  some  rather 
classical  problems  (evaluation  of  integrals)  the  estimated  variance  re¬ 
duction  is  more  easily  found. 

The  evaluation  of  the  Monte  Carlo  method  is  given  below. 

Advantages : 

1)  The  method  is  useful  for  highly  inhomogeneous  media 
and  for  irregular  and/or  higher  dimension  geometries. 

2)  It  is  possible  to  study  perturbations  directly  rather 
than  consider  two  separate  problems. 

3)  By  appropriate  choice  of  the  method  of  analysis  any 
given  property  of  the  shield  can  be  studied,  for 
instance,  reflection  coefficients  rather  than  trans¬ 
mission. 

4)  With  a  proper  selection  of  variance  reduction,  the 
computational  time  may  be  much  smaller  than  direct 
integration  and/or  the  moments  method. 

5)  Very  complicated  interaction  probabilities  are 
readily  incorporated  into  the  collision  mechanics 
without  approximation. 

Disadvantages : 

1)  Frequent  lack  of  reliable  error  estimates. 

2)  Results  may  be  seriously  in  error  without  any 
statistical  indication  available. 
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Disadvantages,  Cont.: 

3)  For  some  problems  the  time  of  computation  may  be 
excessive  -  this  is  usually  true  for  problems  of 
small  asymmetry  and  few  dimensions . 

4)  Distribution  functions  to  be  sampled  may  be  very 
complex  and  time  consuming  to  select  from 

5)  Problem  must  usually  be  designed  to  yield  only 
limited  data  due  to  time  and  storage  problems. 

3 .4  Method  of  Successive  Scatterings 


The  method  of  successive  scatterings  was  developed  for  y-ray 
penetration  problems,  particularly  in  slab  geometry.  The  technique  has 
the  virtue  of  being  applicable  to  multi-layer  shields.  Further,  the 
method  is  only  used  for  finite  shields  which  may  have  any  desired  thickness. 

The  basis  of  the  method  is  to  consider  the  fractional  transmission 
of  photons  of  0,  1,  .  .  .  k  scatterings  within  the  shield.  If  N^(a  ,XQ ,  |i0) 
is  the  fraction  of  the  incident  photons  transmitted  which  undergo  k  collisions 
then 


N(a,X0.^o)  =  4  \(a.X0,M-0) 


where  N  is  the  total  fraction  transmitted.  The  difficult  part  of  the 

method  is  the  determination  of  the  k  >  0.  For  k  *  0,  we  have 

a  a 
o 

N  (a,  X  ,|0r  )  «  e 

O  0  0 

The  factor  N]_  can  be  found  as  follows.  consists  of  all  photons  which 

undergo  a  first  collision  and  then  escape.  The  probability  of  traveling 
a  distance  x  into  the  shield  is  merely 


a  x 


The  probability  of  undergoing  a  collision  within  dx  is  crQ  — .  The 
probability  of  scattering  through  the  angle  0  with  respect  to 

the  initial  direction  is 

as  (9,\o)dO 


cr 
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Finally ,  the  probability  of  being  transmitted  without  further  collision  is 
N  (a  -  x,\,pO  where  X  is  the  wavelength  of  the  scattered  photon  and  ji  the 
direction  cosine.  Thus, 


dNx  =  e 


v 

U  a  dx  a  (0,\  )  dO  N  (a-x 
0  o  s  ’o  -  O  5  5rv 


a  x 
o 


N1(a’Xo^o) 


J  dx  J  dQ  -  as(0,\o)  Nq (a-x,\,n) 


o  Q 


(20) 


The  only  elements  of  the  integral  over  the  solid  angle  which  can  contribute 
to  the  transmission  are  those  for  which  0  such  that  >  0.  We  then  have 

a  X  r  s 

£l£;2Sl 

^ 

-  ct  (Q.X  )  e  (21) 

<  S  ?  O 


The  integration  over  x  is  performed  analytically.  The  resultant  distribution 
is  then  numerically  integrated  over  9;  cp  from  the  known  Klein-Nishina  cross- 
section  „ 


Ni<a>V>V 


dx 

O  t 


dQ 


>  0 


A  completely  analogous  procedure  is  used  for  higher  k.  values. 

Thus  N,  (a,\  ,[Jb  )  is 
k  o  o 


axi 


k~l  a.x. 
-  x  1 


;(a’X0^0>  =  I  I  ^  J  J  d%  Tl  e_  °S(9>X)  * 

o  o  1=0 


o  o 
k  terms 


k  terms 


CTi (a-xi) 


*  e 


“i 


(23) 


In  order  to  make  the  integrations  over  angle  manageable  it  is  assumed  that 

only  forward  scattered  photons  contribute  sensibly  to  the  transmission 0 

Consequently  all  the  \i,  >0.  Although  the  x  integrations  are  analytic,  the 

angular  integrations  are  numeric.  In  particular  2k  integrations  are  needed 

to  evaluate  N.  . 

k 
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Because  of  the  involved  integrations  usually  the  series  are 
terminated  at  k  ^  3,  and  remaining  values  estimated.  Without  considering 
the  details  further,  it  is  evident  that  the  method  has  serious  drawbacks 
for  compact  power  plants.  The  relative  merits  are  listed  below. 

Advantages : 

1)  Useful  for  quick  surveys  for  thin  shields . 

2)  Can  be  modified  to  yield  energy  transmission 
rather  than  number  transmission. 

3)  The  method  has  been  extended  to  multi-layered 
shields  . 

Disadvantages : 

1)  Computational  burden  is  large  for  large  number 
of  scattering  components. 

2)  The  results  do  not  give  all  information  desired, 
i.e.,  the  spectrum  and  direction  of  emergent 
photons  . 

3)  Other  geometries  are  very  difficult. 

4)  The  errors  are  large  for  thick  shields.  In 
particular,  it  is  difficult  to  even  estimate 
the,  errors  . 

5)  Although  the  method  is  actually  an  approximate 
solution  of  the  integral  transport  equation,  the 
approach  is  such  that  much  useful  information  is 
not  available,  in  distinction  to  the  usual  solution 
of  the  transport  equation. 

3.5  The  Removal  Cross  Section  Method 


The  simplest  method  of  treating  neutron  penetration  problems  is 
the  removal  cross  section  method.  The  basis  of  the  method  is  the  observation 
that  a  neutron  collision  with  hydrogen  produces  a  lower  energy  neutron  which 
does  not  penetrate  much  further  due  to  the  increasing  hydrogen  cross-section. 
Thus,  the  first  collision  density  determines  approximately  the  penetration 
properties  of  hydrogenous  shields.  For  a  sufficiently  thick  shield,  the 
flux  density  should  ultimately  become  exponential  in  nature. 

The  method  is  limited  to  neutrons  from  the  fission  spectrum. 

After  a  sufficient  thickness  of  hydrogenous  media  (usually  water)  the 
lower  energy  neutrons  are  thermalized  (and  easily  absorbed)  .  The  very 
high  energy  end  of  the  spectrum  is  attenuated  by  the  fission  spectrum 
itself.  The  result  is  the  spectrum  "hardens"  up  to  a  certain  thickness 
and  then  behaves  roughly  as  a  mono energetic  beam. 
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If  a  slab  of  material  is  introduced  between  the  source  and  the 
hydrogenous  medium,  the  behavior  deep  within  the  shield  (i.e.,  the  water) 
will  be  the  same  as  without  the  slab  but  with  reduced  magnitude.  Let 
cp  (x)  be  the  hardened  flux  distribution  in  the  absence  of  the  slab,  and 
cp^x)  the  distribution  with  the  slab.  For  sufficiently  large  x,  it  is 
phenomenologically  true  that 

cp  (x)  =  acpQ  (x) 

-a  t 
r 

where  O'  may  be  written  e  ,  with  t  the  slab  thickness.  The  coefficient 
ar  is  the  slab  removal  cross-section  for  the  material  and  is  independent 
of  energy,  x,  and  t. 

The  fact  that  such  a  simple  approximation  is  valid  depends 
crucially  upon  having  a  fission  spectrum  and  a  thick  hydrogenous  shield. 

The  measurement  of  cp(x)  (or  90 GO )  is  actually  very  difficult  and  implicit 
in  the  removal  cross-section  measurement  is  the  assumption  that  the  thermal 
neutron  flux  parallels  the  fast  flux  for  sufficiently  large  x. 

The  removal  cross-section  method  can  easily  be  extended  to  cover 
distributed  shield  material  mixed  with  water.  Although  the  method  is  very 
simple  the  actual  physical  measurement  of  ar  for  either  slab  or  distributed 
shield  material  is  rather  difficult. 

From  this  brief  sketch,  it  is  possible  to  rate  the  method. 

Advantages : 

1)  Simple  to  use. 

2)  Can  be  very  accurate  under  appropriate  conditions. 

3)  Correction  for  irregularities  in  geometry,  channels, 
piping,  etc,  are  simple  in  this  model. 

Disadvantages : 


1)  Requires  a  thick  hydrogeneous  shield, 

2)  Lack  of  a  theoretical  model  for  predicting 
removal  cross-section. 

3)  Uncertainties  in  cross-sections  yield  large 
uncertainties  in  penetration. 
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4.0  A  PROPOSED  SHIELD  CALCULATION  PROCEDURE 

For  all  of  the  methods  reviewed  in  Section  3.0,  only  the  numerical 
integration  of  the  transport  equation  gave  all  of  the  desired  information 
regarding  flux  distributions.  The  major  drawback  with  the  method  was  the 
machine  time  required  to  obtain  the  desired  information.  If  the  calcu- 
lational  procedure  could  be  reduced  then  the  method  would  be  the  most 
general  and  most  exact  procedure  available,  save  for  analytic  solutions. 

The  purpose  of  this  section  is  to  outline  a  new  method  of  solution  of 
boundary  value  problems  which  promises  significant  reduction  in  the  time  of 
calculation.  The  method  has  been  applied  successfully  to  diffusion  calcu¬ 
lations.  The  derivation  of  the  technique  will,  therefore,  be  for  diffusion 
calculations.  A  possible  method  of  extension  for  transport  problems  is 
then  discussed.  It  should  be  clearly  understood  that  the  proposed  extension 
has  not  been  attempted  and  hence  constitutes  a  conjectured  procedure.. 

For  illustrative  purposes  the  diffusion  equation  in  the  form 

V*D(x,y)  V<P  +  X  p  (x,y)  cp  =  0  (24) 

is  considered.  The  equation  applies  in  the  rectangular  region  0  £  x  ^  a, 

0  £  y  £  b  .  A  rectangular  network  is  superimposed  upon  the  region  as 
shown  in  Figure  1. 


Figure  1 

The  Difference  Mesh  for  the  Rectangular 


Region  O^x^a;  O^y^b 
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The  diffusion  equation  (24)  is  assumed  to  obey  boundary  conditions  of 
the  form 

9(0, y)  =  f(y), 

cp(x,0)  =  cp(x,b)  =  0  (25) 

Cp(a,y)  =  0 

The  differential  equation  may  be  reduced  to  a  difference  equation  of  the 
form 


a,  i  i  ,  +  b.  .  cp.  -  +  c.  ,  Cp.  +  d.  .  Cp.  .  i  +  e*  ,  <P.  1^.1  55  0  (26) 
j  ,k  T  j*Hl,k  j,kTj,k  j,k  Tj-l,k  j,kTj,k-l  j,kTj,k+l 


Equation  (26)  applies  of  0  <  j  <  J,  0  <  k  <  K,  i.e.,  all  interior  points 
of  the  mesh-  The  boundary  conditions  for  the  finite  difference  equation 
become 


'f’o.k  “  fk 
’j.O-’j.K"0 
*J,k  -  ° 


(27) 


The  usual  process  for  solving  (26)  is  iterative-  An  assumed  distri¬ 
bution,  consistent  with  the  boundary  conditions,  is  iterated  by  the 
algorithm 


b. 

J: 


h 


cpP 

ky' 


j+l,k 


cj  ,k^j-l,k  + 


k-1 


•j.k^ 


p  l 
j  ,fcnJ 


(28) 


The  superscript  p  denotes  the  iteration  index.  Other  iterations  are 
possible.  Characteristically,  iteration  procedures  such  as  (28)  requires 
an  amount  of  time  proportional  to  ( JK)  ^  . 

The  new  method  is  based  upon  the  following  procedure.  At  the  column 
j  =  J-l,  define  ^  as 

Vi,i 

9  T_1  9 


tj-l  = 
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Now  choose  K-l  different  vectors  as  follows 

J  *“  i- 


(29) 


Thus,  the  set  consists  of  the  set  of  K-l  unit  vectors  and  are  complete 

in  the  space  of  dimension  K-l,  The  basic  difference  relation  (26)  is 
factored  in  the  form 


n,l 

5n,2 


n  ,K-1 


cp.  1  ,  =  -  — —  a.  . cp  +  b.  ,cp.  +  d.  ,cp.  +  e.  ,cp.  ..J 

c  .  L  j ,kT j+l,k  J |k  ] ,k  j,kTj,k-l  j,kTj,kfU 

J  5  K 


(30) 


Then  each  vector  is  extended,  by  equation  (30)  ,  through  the  mesh  to  the 

column  j  =  0 ,  Thus ,  K-l  vectors  are  generated  at  column  0,  say  The 

set  of  vectors  is  complete  and  hence  the  boundary  condition  r.  may  be 

expanded  as  ^ 


f . 
J 


(31) 


The  same  expansion  for  the  vectors  then  yields  the  desired  solution. 

The  entire  procedure  requires  (K-l) (j-1)  steps  . 

,  v  There  are  two  objections  to  the  method  as  it  stands.  First,  the  set 
)  may  not  be  complete  and  hence  the  expansion  (31)  is  invalid.  In 
general,  this  is  indeed  the  case.  The  reason  is  that  the  march-out, 
equation  (30)  is  unstable.  That  is.  errors  in  the  vectors  are  ampli¬ 

fied  to  such  an  extent  that  all  Jv)  become  proportional  to  one  another 
irrespective  of  the  value  of  ij_q*  This  follows  since  the  eigenvalues 
of  the  march-out  operator  are  greater  than  unity  in  general. 

However,  steps  may  be  taken  to  prevent  divergence  of  the  march-out. 

In  particular,  if  the  set  ,  which  are  orthogonal  at  j  =  J-l,  are 

periodically  orthogonalized^  then  the  growing  components  are 
filtered  and  the  completeness  of  the  set  is  assured.  The  only  exception 
is  for  X ^  an  eigenvalue  of  the  equation,  in  which  case  no  solution  exists 
in  any  event.  (Incidentally  the  method  can  therefore  be  used  to  solve 
eigenvalue  problems  also.) 


TM-3522 


28 


The  or thogonalization  requires  a  number  of  steps  proportional  to 
(K-l)^  and  the  entire  procedure  requires  (J-l) (K-l)^  steps.  For  particular 
problems  there  may  be  a  significant  reduction,  i.e.?  when  J  >  K.  However, 
the  method  may  still  be  improved. 

The  initial  orthogonal  set  chosen  was  the  unit  vectors.  Alternative 
sets  may  be  used,  for  instance,  the  finite  Fourier  harmonics.  If  the 
solution  is  reasonably  smooth  only  a  few  Fourier  harmonics  are  necessary 
to  specify  the  solution.  Suppose  the  first  M  harmonics  are  sufficient. 

The  number  of  steps,  including  the  re~orthogonalization ,  is  then  (K.-1)  (J-l)  (M) 
The  use  of  an  incomplete  set  of  harmonics  introduces  a  truncation  error. 
However,  for  most  problems  the  amplitudes  of  the  harmonics  decreases  rather 
rapidly  on  either  side  of  the  fundamental.  Of  course,  other  orthogonal 
vectors  than  the  Fourier  series  are  possible  for  use. 

A  practical  procedure  for  selecting  the  number  and  range  of  vectors  is 
to  expand  the  boundary  condition  f^  and  select  only  the  significant  vectors. 
Any  desired  order  truncation  error  may  be  found. 

All  of  the  above  ideas  have  been  successfully  applied  to  the  diffusion 
equation.  It  has  been  found  that  the  method  permits  a  significant  time 
saving  for  solution  of  the  equation.  In  particular,  the  criticality  problem 
has  been  solved  for  the  first  5  critical  eigenvalues  with  ease,  a  problem 
that  is  very  difficult  by  other  means .  The  use  of  a  truncated  orthogonal 
series  has  been  used  to  solve  inhomogeneous  problems  and  experimentally  it 
is  found  that  3  vectors  are  sufficient  to  reduce  the  truncation  error  to 
below  1/2%.  Further  research  is  needed  to  find  more  general  rules  for 
selecting  the  proper  vectors. 

The  application  of  the  method  to  the  transport  equation  is  not  as 
clear  cut  as  the  above.  The  particular  problem  of  concern  is  the  boundary 
conditions.  Consider  a  finite  mesh  in  x-p  spare  as  shown  in  Figure  2. 


Figure  2 

A  Discrete  Mesh  in  r~p  Spare  for  the  1  Group 


Transport  Equation 
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The  appropriate  boundary  conditions  at  x  *  0  include  the  incident 
current,  i.e.,  >  0.  However,  nothing  is  known  about  the  emergent  current, 

i,.e  . ,  u  <  0.  Likewise  at  x  =  a,  the  incident  current  is  zero  for  a  vacuum 
interface.  This  specifies  cp(a,p«),  \i  <  0.  Again,  nothing  is  known  about 
the  flux  for  \i  >  0,  that  is  the  leakage.  The  boundary  condition  at  =  +  1 
for  all  x  are  merely  $2.  =  0,  that  is  no  gradient  of  the  directional  flux.. 

One  possible  approach  is  to  consider  the  flux  density  to  vanish  at  the 
extrapolated  end  point  of  the  region  at  x  *  a...  For  shielding  studies  it  is 
not  clear  that  such  an  approximation  is  proper,  since  the  exit  current  is 
precisely  the  desired  information.  Furthermore,  for  different  energy 
groups,  the  end-point  is  variable.. 

An  alternative  procedure  is  to  divide  the  mesh  at  [i  =  0  into  two 
problems  with  the  condition  that  Cp(x,0)  be  continuous.  The  difficulty  here 
is  the  treatment  of  the  line  at  ^  =  0,.  Further,  the  expansion  vectors  along 
the  portion  of  the  problem  for  >  0  need  not  be  the  same  for  \i  <  0. 

It  is  clear  from  the  above  discussion  that  some  effort  should  be 
expended  in  adopting  the  procedure  to  transport  problems.  The  promise  of 
the  reduction  in  time  justifies  considerable  effort  in  this  direction.. 

Under  the  assumption  that  the  method  can  be  adopted  to  the  transport 
equation,  the  numerical  solution  would  then  be  the  preferred  method  of 
attack  for  compact  power  reactor  shields,.  For  the  very  difficult  geometries 
such  as  rocket  vehicles,  the  numerical  integration  solution  over  the  primary 
shield  must  be  coupled  with  another  method  for  solving  for  scattered  radiation 
into  the  payload.  Monte  Carlo  is  at  present  the  only  reasonable  procedure. 
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